Burgers' equation

In [1]:
%matplotlib inline
In [2]:
%config InlineBackend.figure_format = 'svg'
from utils.snapshot_widgets import interact
from ipywidgets import widgets
from ipywidgets import FloatSlider, fixed
from exact_solvers import burgers
from exact_solvers import burgers_demos

Will create static figures with single value of parameters

In this chapter, we study a simple scalar nonlinear conservation law: Burgers' equation. Burgers' equation models momentum transport in a fluid of uniform density and pressure, and it is the simplest equation that captures some key features of gas dynamics or shallow water flow.

To examine the Python code for this chapter, see:

Burgers' equation has been used extensively for developing both theory and numerical methods, and it will allow us to further explore the Riemann problem for nonlinear conservation laws. The Burgers' equation is given by \begin{align} q_t + \left(\frac{1}{2}q^2\right)_x = 0. \label{burgers0} \end{align} It is a scalar conservation law with the nonlinear flux term $f(q)=q^2/2$. The quasilinear form is obtained by applying the chain rule to differentiate the flux term:
\begin{align*} q_t + qq_x = 0. \end{align*}
This equation looks very similar to the advection equation, with the difference that the advection speed at each point is given by the solution $q$. The nonlinearity in the Burgers flux term is essentially the same as in the momentum conservation equation in fluid dynamics, so studying its dynamics provides a guideline to understand the more complex nonlinear systems arising in the study of shallow water or Euler equations, as we pursue further in Shallow_water and Euler.

In Traffic_flow we consider another scalar nonlinear conservation law that has a similar structure to Burgers' equation and a nice physical interpretation.

Shock formation

As the speed of propagation in Burgers' equation is given by $q$ itself, the peak of a wave travels faster than the rest. This behavior is illustrated in the next figure. The dashed line shows the initial condition while the solid lines show the solution at later times.

In [3]:
interact(burgers_demos.multivalued_solution, 
         t=FloatSlider(min=0.,max=9.,value=7.75), fig=fixed(0));

Notice that at first $q$ remains single-valued for every $x$. However, after some time the crest of the wave overtakes the leading edge. After this time, we obtain a triple-valued solution for certain values of $x$, which is unphysical since momentum is single-valued. The first time this overtaking happens is referred to as the breaking time -- a reference to waves breaking on the beach. It is also the point where the conservation law, in its differential form, breaks down and where the characteristics cross each other for the first time. When characteristics cross, a shock wave forms. Mathematically, imposing a shock will avoid the problem of the solution being multivalued. Where should the shock be located?

If we replace part of the multivalued solution interval with a shock, some mass will be removed (area $A_1$ in the figure below) and some mass will be added (area $A_2$ in the figure below). In the live notebook, the figure below shows possible locations for the shock at a given time; in order to maintain conservation of the integral of $q$, the shock must be placed so that these two areas are equal.

In [4]:
interact(burgers_demos.shock_location, 
         xshock = FloatSlider(min=6,max=10,step=0.25,value=7.75),
         fig=fixed(0));

This geometric reasoning provides a nice intuition for the shock location, but is cumbersome in practice. To determine the path of the shock we can use the Rankine-Hugoniot jump condition, which we will derive in Traffic_flow, and requires that the jump in flux across a propagating shock must be related to the jump in $q$ by $$ s(q_r-q_l) = f(q_r) - f(q_l), $$ at each instant in time, where $s$ is the shock speed at this time.

Since the flux for Burgers' equation is $f(q) = q^2/2$, this gives \begin{align*} s(q_r-q_l) & = \frac{1}{2} (q_r^2 - q_l^2) \\ \Rightarrow \ \ \ \ s &= \frac{1}{2}(q_l + q_r). \end{align*}
In general (as in the image above) the states $q_l$ and $q_r$ just to the left and right of the shock will not be constant and the speed of a shock will change in time.

Shock solution

For Burgers' equation, (or any scalar hyperbolic conservation law with a convex flux function, such as the traffic flow problem considered in Traffic_flow), the solution of the Riemann problem will always consist of a single wave, which may be either a shock or a rarefaction. The analysis above already yields the solution to the Riemann problem in the case that the resulting wave is a shock. The entropy condition, as we will explain below, indicates that a shock will form when $f'(q_l) > f'(q_r)$; i.e. when $q_l>q_r$. In this case, the solution is \begin{align*} q(x,t) = \begin{cases} q_l \quad \text{if} \ \ x/t<s \\ q_r \quad \text{if} \ \ x/t>s. \end{cases} \end{align*}
Below, we plot the solution of the Burgers' equation for an initial condition that leads to a shock (i.e. with $q_l>q_r$).

In [5]:
interact(burgers_demos.shock(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Rarefaction wave

In the previous figure with the hump as the initial condition, we observed that a shock formed on the right side of the hump. However, on the left side, the characteristics spread out and will never cross. This part of the solution is therefore a rarefaction wave. This is the kind of behavior we will observe in the solution of the Riemann problem when $q_l<q_r$.

In the next figure, we consider such a Riemann problem. As time evolves, a rarefaction forms following the advection speed $q$ given by the quasi-linear equation. Therefore, for time $t$, the solution must propagate a distance $x=qt$, so the solution along the rarefaction is then $q = x/t$. As $q_l<q_r$, the smallest and largest displacements are given by $q_l t$ and $q_r t$, respectively.

In [6]:
interact(burgers_demos.rarefaction_figure, 
         t = FloatSlider(min=0,max=9,value=5));

The full rarefaction solution for the Burgers' Riemann problem is then simply given by \begin{align*} q(x,t) = \begin{cases} q_l, \quad \text{for} \ \ x<q_l t \\ x/t, \quad \text{for} \ \ q_l t \le x \le q_r t \\ q_r, \quad \text{for} \ \ x>q_r t. \end{cases} \end{align*}

As we will see in Traffic_flow, the rarefaction solution is always a self-similar solution. This means that it can be expressed as a function of the ratio between position and time $q(x,t) = \tilde{q}(x/t)$, so it remains the same when rescaling both $x$ and $t$ by the same factor. This is a consequence of $q$ being the same along the characteristic in the rarefaction fan. In the Burgers' equation, the form of the rarefaction is particularly simple since the advection speed is simply $q$.

Rarefaction solution

In the figure below, we plot a solution of the Riemann problem with $q_l<q_r$. In this case, we observe a rarefaction.

In [7]:
interact(burgers_demos.rarefaction(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Weak solutions

As we mentioned before, the differential form of the equation breaks down in the presence of shocks/discontinuities. However, the integral form of the conservation law remains valid. We will derive a specific form of the integral conservation law that is particularly useful to prove mathematical results. We first integrate the general conservation law $q_t+f(q)_x=0$ from $x=x_1$ to $x=x_2$ and $t=t_1$ to $t=t_2$ \begin{align} \int_{t_1}^{t_2}\int_{x_1}^{x_2} [q_t+f(q)_x] dx dt = 0, \label{eq:Burgersintclaw} \end{align} This integral can be rewritten in terms of an indicator function $\phi(x,t)$ with compact support in $[x_1,x_2]\times[t_1,t_2]$ (defined to be 1 in this region, zero elsewhere): \begin{align*} \int_{0}^{\infty}\int_{-\infty}^{\infty} [q_t+f(q)_x]\phi(x,t) dx dt = 0. \end{align*}

We can generalize this result by replacing $\phi(x,t)$ by a smooth function with compact support on some region of the $x-t$ plane. In this case, integrating by parts yield \begin{align*} \int_{0}^{\infty}\int_{-\infty}^{\infty} [q\phi_t+f(q)\phi_x] dx dt = -\int_{-\infty}^{\infty}q(x,0)\phi(x,0)dx, \end{align*} where now the derivatives are on $\phi(x,t)$ and not on $q$ or $f(q)$, so the equation still makes sense for discontinuous $q$. The function $q(x,t)$ is called a weak solution of the conservation law if this equation holds for all functions $\phi(x,t)$ that are continuously differentiable and with compact support (bump functions). Note this rules out the original interval chosen in Eq. \ref{eq:Burgersintclaw} since its $\phi$ would not be smooth. However, we can approximate this function arbitrarily well by a smooth function. Note any weak solution is a solution of any form of the integral conservation law and vice versa.

Weak solutions are thus allowed to include discontinuities. The shock solution presented above, for instance, is a weak solution of the Burgers' equation. After characteristics cross, we cannot have strong solutions of the conservation law, and we must resort to weak solutions. However, a potential problem with the notion of weak solutions is that in some cases the weak solution is not unique and it is possible to determine more than one function that satisifes the condition above. This generally happens when discontinuous initial data should lead to a rarefaction wave solution but there is also a discontinous weak solution that is allowed mathematically but is not physically correct.

An unphysical weak solution

We previously mentioned that for Burgers' equation, when $q_l<q_r$, we expect a rarefaction to be the correct solution. This is because if we smeared out the initial data just slightly then there would be smoothly varying characteristics emerging from each point $x$ and these characteristics would spread out. However, for exactly discontinuous initial data, we could also assume the solution consists of a propagating discontinuity and then use the equation in its conservative form to derive the Rankine-Hugoniot conditions, as done in previous chapters and in (LeVeque, 2002). This is mathematically correct and one obtains the same relation for the shock speed as before, $s = (q_l + q_r)/2$ , so the solution is simply \begin{align*} q(x,t) = \begin{cases} q_l \quad \text{if} \ \ x/t<s \\ q_r \quad \text{if} \ \ x/t>s. \end{cases} \end{align*} This unphysical solution is also referred to as an expansion shock, and it is also a weak solution to Burgers' equation. In the next figure, we plot this weak solution.

In [8]:
interact(burgers_demos.unphysical(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Note the speed of the characteristics with respect to the shock. The spreading characteristics clearly indicate that the correct solution should be a rarefaction. In order to be able to specify which of the weak solutions is physically correct, we need to derive a mathematical condition from our physical intuition gained from observing the behavior of the characteristics. This condition is referred as the entropy condition.

The entropy condition

Any condition that allows us to determine whether a weak solution to a conservation law is actually physically correct is called an admissibility condition or more often an entropy condition. This name comes from gas dynamics, where the physical entropy must increase in the gas passing through a shock, according to the second law of thermodynamics. A discontinuity that violates this is non-physical. A mathematical version of an "entropy function" can be defined for other conservation laws. More discussion of this and several other formulations with more detailed explanations are available in the literature, e.g., (LeVeque, 2002).

In the context of the Riemann problem, the entropy condition allows us to determine if the physical solution should involve a shock or a rarefaction. In the case of scalar conservation laws, the entropy condition is quite simple and can be formulated in terms of the flux function of the conservation law as follows: the solution of a scalar Riemann problem will consist of a shock only if $$ f'(q_l) > f'(q_r). $$ In other words, the solution is a shock only if the characteristics are going into the shock as time progresses. If the characteristics are spreading out/going out of the shock as in the last example, one would naturally expect a rarefaction to be the correct solution. In the case of Burgers' equation, the flux function is $f(q)=q^2/2$, so the correct solution is a shock only if $$ q_l > q_r, $$ which can be clearly observed in the interactive solution below. In later chapters, we will see how this condition generalizes to systems of conservation laws.

Convexity

The flux $f(q) = q^2/2$ for Burgers' equation is a convex function, since $f''(q) = 1 > 0$ for all values of $q$. This means that as $q$ varies between $q_l$ and $q_r$ the the characteristic speed $f'(q) = q$ is either monotonically increasing (if $q_l < q_r$) or monotomically decreasing (if $q_l > q_r$). Because of this the solution is always either a single rarefaction wave or a single shock wave.

The flux $f(\rho) = \rho(1-\rho)$ simple traffic flow model considered in Traffic_flow also has the property that $f'(\rho)$ varies monotonically with $\rho$ (in the context of hyperbolic PDEs, a flux function is referred to as convex if either $f''(q)\ge0$ for all $q$ or $f''(q)\le 0$ for all $q$). Since $f''$ is always negativefor the traffic flow model, the correct Riemann solution is a shock if $\rho_l < \rho_r$, or a rarefaction wave if $\rho_l > \rho_r$.

The solution to the Riemann problem can be much more complicated if $f'(q)$ is not monotonically varying between $q_l$ and $q_r$, i.e. if $f''(q)$ changes sign. In this case the Riemann solution can consist of multiple shock and rarefaction waves. The nonconvex case is explored further in Nonconvex_scalar.

Interactive solution and examples

The live notebook next shows a full interactive solution for the Riemann problem for Burgers' equation. The values of the initial conditions and the time can be modified to observe their effect on the characteristic structure and the solution.

In [9]:
burgers.riemann_solution_interact()

Animations of more complex solutions

We can now explore a few more examples that are representative of phenomena we will observe in more complicated systems, with more interesting initial data that the single jump discontinuity of a Riemann problem. For the animations shown below, the solution is computed numerically using PyClaw.

Example 1: Bump initial condition

The animation below, or on this webpage, shows the solution of Burgers' equation for an initial Gaussian hump.

In [10]:
burgers_demos.bump_animation(numframes = 50)
Out[10]:


Once Loop Reflect

Example 3: Three state solutions

When three piecewise constant values are given as initial condition, one can solve a Riemann problem across each discontinuity. For instance, consider the Burgers' equation with its initial condition given by

\begin{align*} q_0(x) = \begin{cases} q_l \quad \text{if} \ \ x < -1 \\ q_m \quad \text{if} \ \ -1\le x \le 1 \\ q_r \quad \text{if} \ \ x > 1. \end{cases} \end{align*}

We can decompose this into two Riemann problems, one at $x=-1$ and another one at $x=1$. Each of these two Riemann problems will yield either a shock or a rarefaction. The interesting part is what happens when the generated shocks or rarefactions interact. Following the animation below, let us assume that both Riemann problems produce a right-propagating shock. However, the shock generated at $x=-1$ propagates faster, so it will eventually reach the shock originally generated at $x=1$. At the point in time when the shocks collide, one can simply restate the problem again as a Riemann problem and solve the whole problem analytically. We again use PyClaw to solve the problem and show an animated solution. The first plot shows the solution $q(x)$ evolving through time, while the second one shows the characteristic structure in the $x/t$ plane, where the shocks are shown as wide lines and the characteristics as thin lines. In this figure, one can clearly see the two shocks collide to become a new one. In the second plot, the time is marked by an horizontal dashed line. Note one can also solve this problem analytically. This animation can be viewed in the live notebook on or this webpage.

In [11]:
burgers_demos.triplestate_animation(ql = 4.0, qm = 2.0, 
                                    qr = 0.0, numframes = 50)
Out[11]:


Once Loop Reflect

As one would expect, this can become more complicated when one or both of the waves generated are given by rarefactions. The animation below, or on this webpage, shows how a right-propagating shock can overtake a rarefaction. However, there are consequences for the shock after interacting with the rarefaction; its speed is changed. This is clear in the $x/t$ plane, where the slope of the shock is changed after interacting with the rarefaction. Note the shock is again marked as a wider line.

In [12]:
burgers_demos.triplestate_animation(ql = 4.0, qm = -1.5, 
                                    qr = 0.5, numframes = 50)
Out[12]:


Once Loop Reflect

Analogously, one can also have a rarefaction colliding into a shock. As we can see in the animation below, or this webpage, this can even change the direction in which the shock is propagating. This happens because the velocity of the shock depends on the slopes of the impinging characteristics. When a shock interacts with a rarefaction, the varying slopes of the characteristic curves produce a constant change of the shock propagation velocity. Try modifying the values of $q_l$, $q_m$ and $q_r$ in order to see what other behavior can be observed. Can one observe two interacting rarefactions?

In [13]:
burgers_demos.triplestate_animation(ql = -1, qm = 3.0, 
                                    qr = -2, numframes = 50)
Out[13]:


Once Loop Reflect